#!/usr/bin/perl

# Utilizza il programma HBPLUS che calcola i legami ad idrogeno presenti all'interno di una proteina.
# Il file di output generato riporta:

# 1° colonna: numero residuo donatore
# 2° colonna: atomo donatore
# 3° colonna: numero residuo accettore
# 4° colonna: atomo residuo accettore
# 5° colonna: distanza atomo donatore-atomo accettore in amstong
# 6° colonna: angolo formato tra donatore e accettore con vertice sull'H in gradi

# lo script vuole come input un file .pdb: es ./pass_H_bond.pl 1CTF.pdb
# il file pdb viene passato al programma pdbshifter di VICTOR che rinumera i residui partendo da 1
# viene creato un file pdb temporaneo che alla fine viene cancellato
# tale file pdb viene passato al proframma HBPLUS che produce un file di output contenente tutti i legami ad H tra donatore e accettore e viceversa
# il file viene nominato con estensione .hb2
# tale file viene passerizzato e dalla lista vengono estratti i valori di interesse prima menzionati
# i valori vengono esportati nel file OUTPUT_pass_H_bond
# il file .pdb temporaneo viene poi cancellato
# il file .hb2 viene cancellato 

$hbplus_prg = "./hbplus/hbplus";


($ARGV[0])||(die "Specificare il file PDB\n\n");

qx{pdbshifter --renum -i $ARGV[0] -o TMP_MICHELE.pdb};


print qx{$hbplus_prg TMP_MICHELE.pdb};

open(F, "./TMP_MICHELE.hb2");
open(OUT, ">./OUTPUT_pass_H_bond");
while(my $riga = <F>)
{
    #  A0062-GLY N   -0001-SO4 O1  2.75 MH  -2 -1.00 128.8  2.01 108.3 121.0    24
    #   ----     -     ---     --  ----              -----                        

    if($riga =~ /^[\-A]\d{4}\-\S+\s+(\S+)\s+[\-A]\d{4}\-\S+\s+(\S+)/)
    {
        my ($aa1, $atom1, $aa2, $atom2, $distanza, $angolo) = $riga =~ /^[\-A](\d{4})\-\S+\s+(\S+)\s+[\-A](\d{4})\-\S+\s+(\S+)\s+([\d\.]+)\s+\S+\s+\S+\s+\S+\s+(\S+)/;
        $aa1 =~ s/^0+//gi;
        $aa2 =~ s/^0+//gi;
        #print "$riga";
        print OUT "$aa1\t$atom1\t$aa2\t$atom2\t$distanza\t$angolo\n";
	
    }

}
close F;
close OUT;
qx{rm TMP_MICHELE.pdb}; 
qx{rm TMP_MICHELE.hb2};

